SYDE 556/750: Simulating Neurobiological Systems

Terry Stewart

Accompanying readings: Chapters 4, 5

Temporal Representation

• In the previous lecture, we covered basic representation
  • take a vector $x$ and convert it into a bunch of neural activity values $a_i$, and then convert back to the original value $\hat{x}$.
• What happens if $x$ changes over time?
  • i.e. if we have $x(t)$, how do we get $a_i(t)$?
  • seems pretty easy:
    • $a_i=G[\alpha e \cdot x + J^{bias}]$
    • $a_i(t)=G[\alpha e \cdot x(t) + J^{bias}]$
    • where $G[J]= {1 \over {\tau_{ref}-\tau_{RC}ln(1-{1 \over J})}}$

• Unfortunately, that formula for $G[J]$ is an average firing rate over a long period of time.
  • What does it mean to be firing at 10Hz for 10ms?
• We need to model what's actually happening to that neuron over time

Leaky Integrate-and-Fire neurons

• Pretty much the standard neuron model
  • simple
  • produces spikes
  • limiting case of other more complex neuron models
  • all parameters map onto known properties of real neurons

• bilipid cell membrane acts as a capacitor
• flow of ions through the cell membrane ion channels is a resistor
  • this is the "leak" current
• resistance and capacitance together give $\tau_{RC}$
  • $\tau_{RC}=R C$
• when voltage passes some threshold, neuron emits a "spike"
  • very fast, sharp response that's pretty much the same each time
  • spike at time $t_n$ modelled as $\delta(t-t_n)$

• neuron recovers over a short period of time

  • refractory time $\tau_{ref}$
  • during this time the voltage is reset to some resting level

• typical values

  • $\tau_{RC}$: 0.02 seconds
  • $\tau_{ref}$: 0.002 seconds
  • reset voltage (a.k.a. resting potential): -70mV
  • firing threshold (a.k.a. threshold potential): -55mV
    • NOTE: we normalize this so that reset is at 0 and firing is at 1
    • this does not affect anything about the behaviour of the model

• For lots more info, see Unit 2 in this online course

The RC circuit

• let's ignore the spike for now

• Kirchhoff's Current Law
  • Currents at a point have to add up to zeros
  • Current coming into the neuron has two paths
  • $J_{input} = J_R+J_C$
• Ohm's Law
  • Formula for current through a resistor
  • $ J_R = {V \over R}$
• Capacitor
  • $ C = {Q \over V}$
  • $ Q = C V $
  • $ {{dQ} \over {dt}} = C {{dV} \over {dt}} $
  • $ J_C = C {{dV} \over {dt}} $

$ J_{input} = {V \over R} + C {{dV} \over {dt}} $

$ {{dV} \over {dt}} = {1 \over C}[J_{input} - {V \over R}] $

$ {{dV} \over {dt}} = {1 \over {RC}}[R J_{input} - {V}] $

$ {{dV} \over {dt}} = {1 \over {\tau_{RC}}} [R J_{input} - {V}] $

• What does this look like?

• If $J_{input}$ is constant (and we'll write it as $J$ for simplicity)

  • $ {{dV} \over {dt}} = {1 \over {\tau_{RC}}} [R J - {V}] $

  • $ V(t) = JR(1-e^{-t/\tau_{RC}}) $

• So how long does it take to get to threshold?

  • $ V_{th} = JR(1-e^{-t_{th}/\tau_{RC}}) $
  • $ t_{th} = -\tau_{RC} ln (1-{V_{th} \over {JR}}) $

• So how long between spikes?

  • $t_{th} + \tau_{ref}$

• Firing rate

  • $a = {1 \over {t_{th} + \tau_{ref}}}$
  • $a = {1 \over {\tau_{ref}-\tau_{RC}ln(1-{V_{th} \over {JR}})}} $

• Simplify by assuming $V_{th}=1$ and $R=1$

  • $a = {1 \over {\tau_{ref}-\tau_{RC}ln(1-{1 \over J})}} $

Spikes

• But we don't want constant $J$.
• We want it to change over time (since $x(t)$ also changes over time)
• So let's do a computer simulation of the differential equation
  • $ {{dV} \over {dt}} = {1 \over {\tau_{RC}}} [R J_{input} - {V}] $

• Let's try a sine wave input

• So this is how we're going to model neurons over time
• Weaknesses of this approach
  • These are "point" neurons (i.e. we're assuming the body of the neuron is just one infinitessimal point, rather than having lots and lots of compartments, each with its own RC circuit)
  • assumes R is a constant
  • assumes $V_{th}$ is a constant
  • assumes no adaptation (it's harder for a neuron to fire after it has recently fired)
  • doesn't consider nonlinearities at the input
• Each of these could be added
  • Adding these increases computational cost
  • Not going to worry about it now
  • But we should make sure everything we do would also work with more detailed neuron models

The time vs. rate-coding debate

• There is a standard ongoing debate in the neuroscience literature about whether neurons use a "rate" code or a "timing" code
• rate coding
  • The important thing is the firing rate over some window of time (~100ms)
• time coding
  • The precise pattern of spike generation carries information
• We don't find this to be a useful distinction
  • Some people call what we do rate coding, some call it time coding
• The important thing is how to do the decoding...

Temporal Decoding

• So that's what happens when we allow $x(t)$ to vary over time
  • We keep $J(t)=\alpha e \cdot x(t) + J^{bias}$
  • Feed that into the neuron model $G[J]$
  • We get out a sequence of spikes

• How do decode this?

  • How do we get an estimate of $x$ given that sequence of spikes?

• Can we stick with linear decoding?

  • Let's consider a case with 2 neurons
  • For simplicity, use the same $\alpha$ and $J^{bias}$, but use $e_1=1$ and $e_2=-1$

• Can we still decode $x$ using a sum of activities?

  • $\hat{x}(t)=\sum_i{a_i(t) d_i}$

• Let's use exactly the same technique we did before to find $d$

  • $ d = \Gamma^{-1} \Upsilon $
  • $ \Upsilon_i = \sum_x a_i x dx$
  • $ \Gamma_{ij} = \sum_x a_i a_j dx $

• What if we add more neurons?

• What's happening here?
• At a given instant in time $t$, only a very small number of neurons are firing
  • if no neurons are firing at that instant, $\hat{x}=0$
  • there's usually only one or two neurons firing at any given time
  • as $dt$ gets smaller, the problem gets even worse
• What can we do?

Temporal filter

• We want a spike to contribute to other points in time
• So we want to convert a train of spikes into some continuous function

• Convolution:

  • $\hat{x}(t)=\sum_i{a_i(t) * h(t) d_i}$

• What should we use for $h(t)$?

  • Maybe a gaussian?
  • $h(t) = e^{-t^2 / {2\sigma^2}}$

• Better, but how can we improve this?
  • Different shapes for $h(t)$?
  • Different $\sigma$?
• Is there a general solution?

Optimal filters

• Let's try to find the optimal filter for a particular special case
  • Two neurons
  • same $\alpha$ and $J^{bias}$ but different $e$ (-1 and +1)

$\hat{x}(t)=\sum_i{a_i(t) * h(t) d_i}$

• Since there are two neurons and they are opposites of each other, we can assume that $d_1=-d_2$ (i.e. they both contribute equally, but oppositely)

  • $\hat{x}(t)= \sum_i{a_i(t) * h(t) d_i} $
  • $\hat{x}(t)= a_1(t)*h(t)d_1 + a_2(t)*h(t)d_2$
  • $\hat{x}(t) = (a_1(t)-a_2(t))*h(t) d_1$
  • $\hat{x}(t) = (a_1(t)-a_2(t))*h(t)$ (since we can roll the constant into $h(t)$

• Let's call $a_1(t)-a_2(t)=r(t)$, for response

$\hat{x}(t)=r(t) * h(t)$

Random inputs

• Can't just optimize for decoding a sine wave

  • needs to work for any random input

• What do we mean by a random input?

• Maybe that's a bit too random
• Real signals don't change that quickly
  • how do we get a random signal that doesn't change quickly?
  • i.e. a signal with a limit on its maximum frequency?

• Generate random Fourier coefficients and convert back to time domain
  • randomly generate coefficients from a Gaussian distribution
  • for $\omega$ values outside of some range, set to zero

• Things to remember
  • $X(\omega)$ values are complex numbers
  • When we plot $X(\omega)$ we actually plot $|X(\omega)|$
  • Since $x(t)$ should be real (i.e. no imaginary component), we need to constrain $X(\omega)$
  • This constraint is that $X(\omega)=X(-\omega)^*$
  • (the complex conjugate: the real parts are equal, and the imaginary parts switch sign)

Finding an optimal decoder

• $\hat{x}(t)=r(t) * h(t)$
  • $r(t)$ is the response of the two neurons combined together
  • $r(t)=a_1(t)-a_2(t)$
  • $h(t)$ is the filter we are trying to find
• how do we find $h(t)$?
  • convolution is annoying; let's get rid of it
  • one nice thing about convolution is that it turns into multiplication when you do a Fourier transform

$\hat{X}(\omega) = R(\omega)H(\omega)$

• we want to find $H(\omega)$ that minimizes the error

  • $E = (X(\omega)-\hat{X}(\omega))^2$ (added up over all time)
  • $E = (X(\omega)-R(\omega)H(\omega))^2$ (added up over all time)

• but, we know that our signal $x(t)$ can be written in the frequency domain as a sum of sine waves (that's the Fourier transform)

  • we've already computed those for our signal $X$; that's how we made it in the first place

  • so, we can write our error in terms of different frequency values $\omega$

  • $E(\omega) = (X(\omega)-R(\omega)H(\omega))^2$

  • now we can take the derivative and set it equal to zero to do the minimization

  • $H(\omega)= {{X(\omega)R^*(\omega)} \over {|R(\omega)|^2}} $

  • note the complex conjugate that appears there

• so now we can find $H(\omega)$ given the Fourier transform of our signal $X(\omega)$ and the Fourier transform of the spiking response $R(\omega)$

Notes for comparison to the textbook

• Here we're using the convention that lower-case letters are for time-domain functions and upper-case letters are used for Fourier transfomr
  • $x(t) \rightarrow X(\omega)$
  • $r(t) \rightarrow R(\omega)$
  • $h(t) \rightarrow H(\omega)$
• the textbook uses $A(\omega)$ instead of $X(\omega)$
  • $A$ for amplitude
  • but we're already using $A$ for the matrix for activities of neurons so let's avoid re-using the letter

• Now we can find our optimal filter

  $H(\omega)= {{X(\omega)R^*(\omega)} \over {|R(\omega)|^2}} $

Testing how well it generalizes

• if it's a good decoder, it should work for other inputs too

• It's 5x worse on a different signal
• Not bad, but could use some improvement

Improving the optimal decoder

• We're only training on a pretty small set of data
  • So let's increase $T$
  • This helps, but can get computationally expensive to generate more and more data
• There's also a bit of a shortcut
  • We've already created a pretty long signal ($T=1.0$ seconds)
  • We know the filter should probably decay towards zero
    • A spike shouldn't affect values far away from when it happens
  • So let's take our one existing signal and chop it up into lots of little signals

• Add up a whole bunch of these with different slices
• But, this introduces some weird edge effects at the beginning and ending of the window
• So we should do some sort of smoothing of the data
  • Let's go with a Gaussian

• So we're going to take our original data ($x(t)$ and $r(t)$) and shift it around, multiply it by some other function (the Gaussian), and then add up the results
• There's a word for that
• This is actually a pretty common mathematical operation that we've already used

$H(\omega)= {{(X(\omega)R^*(\omega)) * W(\omega)} \over {|R(\omega)|^2*W(\omega)}} $

How well does this work?

• As expected, it's a bit worse on the original signal
• But it doesn't get hugely worse when tested on a different signal
• What about higher frequency inputs? Or lower frequency?

• the decoder needs to be made for the right range of frequencies

Optimal filters reconsidered

• So we have now found the optimal filter to apply to the spike data to decode out information
• Note a limitation: we've only done this for the case of 2 neurons!
  • But we could use the same filter for many more neurons
• What does this filter look like?

• What is this saying?
• Every time a spike occurs, when we interpret the data we should replace the spike with this shape
• What does this mean?

• If we're analyzing the data after the fact, this is fine
• But what if we're just sitting here obseving the neural activity and wanting to know what's being represented right now?
  • Notice that the blue line (the filtered value we're using for decoding $\hat(x)$) starts going up before the spike happens
  • The filtered value is using information from the future.
• This means that when we look at $\hat{x}$ this way, we're getting a better estimate of $x$ than it is possible for other neurons to get
  • (assuming neurons aren't psychic)
• What could we do instead?

Biologically plausible filter

• We're looking for something that takes a spike and turns it into some smooth signal
• Instead of coming up with the optimal filter, is there anything from neuroscience we could use instead?
• What actually happens when a spike occurs?
  • What is the input to the next neuron that occurs when a spike happens?

• A spike causes a vesicle to release neurotransmitter
• The amount of neurotransmitter causes current to flow into the next neuron
• The neurotransmitter is slowly reabsorbed
• What is the current that goes into the next neuron thanks to a single spike?

(data stolen from here)

• So other neurons don't get spikes as inputs
  • they get these post-synaptic currents

$$ h(t) = \begin{cases} e^{-t/\tau} &\mbox{if } t > 0 \\ 0 &\mbox{otherwise} \end{cases} $$

• this is the standard simple model

  • $\tau$ is different for different neurotransmitters and different neurotransmitter receptors
  • GABA-A-R: $\tau=10.41 \pm 6.16$ ms
  • AMPA-type GluR: $\tau=2.2 \pm 0.2$ ms
  • NMDA-type GluR: $\tau=146 \pm 9.1$ ms

• for more complex models, some people use

$$ h(t) = \begin{cases} t^n e^{-t/\tau} &\mbox{if } t > 0 \\ 0 &\mbox{otherwise} \end{cases} $$

• What about with more neurons?

• what about different neurotransmitter time constants $\tau$?

A shortcut

• What happens if we find our decoders ignoring spikes?
  • Just compute decoders using the non-spiking version of the neurons
  • Then run the model with spikes, filter the data with $h(t)$ and decode

• Why does this work?
• Consider a single neuron with a constant $x$
  • when we found decoders using a rate mode, the firing rate was $a$
  • when we use spikes, we get $a(t)$, which is a sequence of $0$s and $1$s
  • if the rate neuron is a good model of the spiking neuron, then the firing rate (the sum of $a(t)$ divided by the total time) should approach $a$
• but what about $h(t)$?
  • it should just blur $a(t)$, pushing it towards this average
  • but this will only happen if $\int_{-\infty}^{\infty} h(t)dt=1$
  • otherwise $h(t)$ will scale the value as well
• so this trick only works if we normalize $h(t)$

• it acts like an infinite amount of training with constant inputs

  • with more complex neuron models, there may be "memory" effects where the past history of the neuron doesn't end up averaging out like this
  • for those neuron models, there is no rate approximation anyway

• the filtered spiking neuron model gives some sort of variability around the ideal $a$ value

  • this can be treated as noise
  • we can find the variance of this noise (see Appendix C.1)
  • but we're already making sure decoders are robust to noise
  • so we don't have to do anything (except maybe slightly increase the amount of noise the decoders are supposed to be robust to)